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Abstract 

Scatterings of electrons at quasiparticles or photons are very important for many topics in solid 
state physics, e.g., spintronics, magnonics or photonics, and therefore a correct numerical treatment 
of these scatterings is very important. For a quantum-mechanical description of these scatterings 
Fermi’s golden rule is used in order to calculate the transition rate from an initial state to a final 
state in a first-order time-dependent perturbation theory. One can calculate the total transition 
rate from all initial states to all final states with Boltzmann rate equations involving Brillouin zone 
integrations. The numerical treatment of these integrations on a finite grid is often done via a 
replacement of the Dirac delta distribution by a Gaussian. The Dirac delta distribution appears 
in Fermi’s golden rule where it describes the energy conservation among the interacting particles. 
Since the Dirac delta distribution is a not a function it is not clear from a mathematical point of 
view that this procedure is justified. We show with physical and mathematical arguments that this 
numerical procedure is in general correct, and we comment on critical points. 

PACS numbers: 02.60.Cb, 02.60.Jh, 02.60.Nm 
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I. INTRODUCTION 


In solid state physics scatterings of electrons at periodic perturbations (quasiparticles or 
photons) are very important for many research helds and we give three examples in the 
following: 

1. In all-optical switching experiment^ a thin ferrimagnetic him, e.g., GdFeCo, is ir¬ 
radiated by a femtosecond laser pulse which can be linearly or circularly polarized 
and thereafter a demagnetization with subsequent switching of the magnetization can 
be observed under certain preconditions. The fundamental mechanisms are strongly 
debated at the moment, however, electron-photon scatterings, electron-phonon scatter¬ 
ings and electron-magnon scatterings certainly play a big role for the demagnetization 
of the ferrimagnetic him. 

2. In ultrafast demagnetization experimentd^l a thin ferromagnetic him, e.g., Ni or Fe, 
is irradiated by a femtosecond laser pulse which is normally linearly polarized and 
thereafter an ultrafast demagnetization (on the time scale of about 100 femtoseconds) 
without switching of the magnetization can be observed. The magnetization recovers 
on a time scale of several picoseconds. Despite many years of research the fundamental 
mechanisms are still unclear but scatterings of electrons at phonon^^El or at magnon^ 
or at electrond^ have been discussed intensively. 

3. Spin-polarized currents are important for devices in spintronic^, e.g., spin-transistors 
or spin-diodes. The lifetime of the spin-polarized electrons is crucial for the spintronics 
devices. The lifetimes are determined by scatterings of electrons at quasiparticles and 
at interfaces or defects. 

A correct numerical calculation of the various scattering processes is important for the 
understanding of these ehects in solid state physics. In quantum mechanics Fermi’s golden 
rule gives the transition rate j,j^, from an initial electronic state in a solid with energy 
to a hnal electronic state Tyk' with energy e^/k' (jjb band indices; k, k'; wavevectors) 
due to a periodic perturbation arising from a (quasi)particl^ 

l"^jkj'k'l ■ (^i'k' ~ (^jk i hc^q^)) . (1) 
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±hucix is the energy of the involved (quasi)particle (q: wavevector, A: polarization) which 
may be, e.g., photons, phonons, magnons, plasmons etc. with frequency cjqA for absorption 
(plus sign) or emission (minus sign), and scattering matrix element 

where |F) and \F') are the initial and hnal (quasi)particle states and hhqA is the scattering 
operator. Thereby, momentum conservation k ± q = k' + G is demanded (G: reciprocal 
lattice vector). Fermi’s golden rule is the hrst-order approximation of the time-dependent 
quantum-mechanical perturbation theory. It implies that the scattering processes are Marko¬ 
vian which means that a scattering process does not depend on preceding scattering pro¬ 
cesses. Fermi’s golden rule is only valid in a time window where the perturbation time on 
the one hand must be short enough because of the hrst-order approximation and on the 
other hand must be long enough in order to replace the sin(a;)/x-function appearing in the 
derivation of Fermi’s golden rule by the Dirac delta distribution. The validity of Fermi’s 
golden rule for a magnetization dynamics on the 100 fs timescale is critically discussed in 
Ref.H 

Normally, one is not interested in a specihc transition rate ITjkj'k' an initial state 
Tjk to a hnal state but in the total transition rate hFtotai from all initial states to 
all hnal states. Thereby k and k' are related via k ± q = k' -|- G if the scattering is at a 
quasiparticle with wavevector q. This is calculated with Boltzmann rate equationd^ 

VFtotai = (3) 


where 


fkin 




[1 Ujk] fkyk'jk 


IRout 



^ [ d^k [ d^k' [1 - n^k'] 


( 4 ) 

( 5 ) 


Dbz is the Brillouin zone (BZ) volume and n is the distribution function for the electrons. 

Often one is also interested in the rate of change of the distribution function n^k due to 
scattering which is also calculated with Boltzmann rate equation^ 

dujk _ 1 
dt Dbz 


/BZ 


d /c j Uj/k' [1 Ujk] hFj'k'jk [1 /^j'k'] hFjkj'k' 


( 6 ) 
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So we have to calculate Brillouin zone integrals of the form 


/ d3fc9(k)i(e(k)). (7) 

Jbz 

Because the quantities ejk, e^k', hhjvk'jk’ ^jkj'k' t)e calculated numerically only for a 
hnite number of k-points, hnite k-point grids have to be used for the numerical calculation of 
the total transition rate fhtotai or of the rate of change of the distribution function dn^k/dt. 
Thereby, energy conservation e^/k' = £jk i hw^^x and momentum conservation k ± q = 
k' + G have to be fulhlled, however, energy conservation in combination with momentum 
conservation is in general never fulhlled for a hnite k-point grid. Therefore, the Dirac delta 
distribution has to be replaced by a “smeared” delta function in order to obtain a result 
which approximates the integral (which is done, e.g., in Refsand in very many other 
papers) . To do this, often the following equation is used 


d^k g{k) (5(£(k)) 


'BZ 


/ g{k) 

'bz a/TTCT 


exp 


£2(k) 




( 8 ) 


and the smearing parameter a has to be chosen appropriately, see Sec. Ill This means that 
the contribution of a certain grid point to the total transition rate hhtotai or to the rate of 
change of the distribution function drij^/dt is small if the energy conservation is fulhlled 
very badly, and vice versa the contribution is large if the energy conservation is fulhlled very 
well. However, from a mathematical point of view it is not obvious that Eq. ([^ holds since 
the Dirac delta distribution is not a function and the smearing is with respect to the energy 
£ but the integration is with respect to the wavevector k. The problem is explained in more 
detail in Sec. m 

Mathematical proofs of Eq. ([^ under certain preconditions can be found in Ref. flTl 
theorem 7.2.1, and in Ref. [T5l theorem 6.1.5, however, the proofs are for general distributions 
and are very abstract. We want to show in this article that Eq. ([^ is correct by using also 
physical arguments. 

In Sec. we explain in detail the problem which arises when in a Brillouin-zone inte¬ 
gration Dirac’s delta distribution is approximated by a Gaussian. This is done in many 
papers without giving any justification. We therefore think that the outline of this prob¬ 
lem is a novelty per se. Then we give a justihcation of the Gaussian smearing method by 
mathematical and physical arguments. Each of these arguments has been used in other 
contexts in previous papers. The novelty of our paper is that we use these arguments to 
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justify the Gaussian smearing method for the Brillouin-zone integration. First, we consider 
a coordinate transformation from the wavevector variables k = (fcj., ky, k^) to the variables 
e, 'd, ip where e is the energy and'd, ip are variables for the surface of constant energy. This 
transformation involves the Jacobian J {k^, ky, kz)- The inverse function theoremP^says that 
if this Jacobian is nonzero at a k-point, this transformation is invertible. Then the integra¬ 
tion in k-space including 6{e) can be represented by an integration over e, 'd,p oi a. function 
which involves the Jacobian J = [J {k^, ky, kz)]~^ and, which now can without any 

problem be replaced by a Gaussian. The problem is that there are special k-points where 
VkC (k) = 0. For these special points the Jacobian J [k^, ky, kz) is zero, and the reverse 
transformation involving J (e, 'd, p) is not dehned in a rigorous mathematical interpretation. 
According to a general theorem of M. Morse the dispersion relation e (k) exhibits such spe¬ 
cial points because it is periodic in all components. There are special points which can be 
identihed easily, e.g., the F-point and points on the Brillouin zone boundary. These points 
can be avoided by shifting the grid of k-points considered in the Brillouin zone integration 
accordingljff^. Other special points cannot be easily found, and they might be in the shifted 
k-point grid. Van Hove has shownl^ that for three dimensions the appearance of these spe¬ 
cial points does not appreciably modify the result of a numerical integration. We motivate 
these steps by physical reasoning. 


In Sec. Ill we give practical hints for the appropriate choice of the smearing parameter 


a. Finally, our results are summarized in Sec. IV 


II. NUMERICAL INTEGRATION OF THE DIRAC DELTA DISTRIBUTION 


It is very well known that in integrals involving the Dirac delta distribution the distribu¬ 
tion can be replaced by a Gaussian for the limes a —?■ 0. It reads 


+O0 r-\-oo 

d£ g{e) J(£) = lim / de g{e) 


'-0O 


exp — 


(9) 


a/Fct \ a 

where g{e) is a continuously differentiable function which depends on the energy £. The 
Dirac delta distribution is approximated by a Gaussian 


b+OO 


r'+oo 


de g{e) 6{e) ^ / de g{e) ^^exp ( 

TTU \ 


( 10 ) 


for a numerical calculation of the integral and a has to be chosen appropriately, see Sec. Ill 
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However, the integrals in Eqs. ([^, Q, ([^ and (|^ are not over the energy e but over the 
wavevector k. For the sake of simplicity we discuss the following integral 


d^k (7(k) 5(£(k)) 


( 11 ) 


'BZ 


where ^'(k) is a continuously differentiable function of the wavevector k, and the generaliza¬ 
tion to the expression 6 ± kwq\)) used in Eqs. (|^, Q, (|^ and (|^ is straight¬ 

forward. In explicit numerical calculations it is always assumed that also the relation 


d'^k 5 f(k) 5(e(k)) = lim / d^fc ^ffk) ^ exp ^ ^ 

(T^O ' - 


( 12 ) 


JBZ \ 

holds without giving any justihcation, reference or comment and that this may be approxi¬ 
mated by 


d^k g(k) 5(£(k)) 


'BZ 


[ g{k) —^ exp 
Ibz V 


(13) 


However, the Dirac delta distribution is dehned by Eq. 0 and not by Eq. (p(^. We show 


in the following how the use of Eq. (13) can be justihed. 


We consider a coordinate transformation from the wavevector variables k = [kx, ky, kz) 
to the variables e, 'd, ip {e: energy; "d, ip: variables for the surface of constant energy) 


£ = e{kx,ky,kz) 
d i}(^kxj ky^ kz^ 
ip ipij^xi kyi kz) 


(14) 


where the energy dispersion relation e{kx, ky, kz) is known and the surfaces of constant energy 
can be parametrized with two variables d and ip. The inverse function theorem say^ that 
every continuously differentiable, vector-valued function which maps values from an open 
set of M” to other values of an open set of M” (so-called coordinate transformation, e.g.. 


Eq. (14)) and whose Jacobian determinant 


J{kx,ky,kz) = det 


^ de 

de 

de 

dkx 

dky 

dk^ 

&■& 

(9i? 

d'S 

dkx 

dky 

dkz 

d^p 

dip 

dip 

\ dkx 

dky 

dkz / 


(15) 
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is non-zero at a point is invertible in the neighborhood of this point, i.e., the reverse trans¬ 
formation of Eq. (14) 


kx = kx{e,^,ip) 
ky = ky{e,^,(p) 

k^ = k^{e,^,ip) (16) 


exists and can in principle be given in the neighborhood of every point {kx, ky, k^) if the 
above-mentioned conditions are fulhlled. 

If Eq. ( |l4| ) is invertible in the neighborhood of every point k = {kx, ky, kz) —whereby only 
points k with £(k) = 0 are relevant because of the Dirac delta distribution in Eq. ( [IT| — , it is 
possible to make for this neighborhood a local coordinate transformation (using Eq. ([I6|) for 
the function 5 f(k) = g{e, "d, (p) appearing in Eq. (In}. Then, the integral over the wavevector 
k can be replaced by the integral over the variables s,'&,(p 


d^k g{k) h(e(k)) 


'BZ 


= y j d'd j dip\J{e,T9,ip)\-g{kx{e,i9,ip),ky{e,T9,ip),kz{e,T9,ip)) ■ 6{e) 

= [ de [ di^ [ dip\J{e,i^,ip)\-g{e,i^,ip)-5{e) (17) 


where J{e,'d,ip) is the Jacobian determinant of the reverse transformation (16) 


J{e, d, ip) = det 


/ 

dkx 

dkx 

dkx ^ 


de 

d'& 

dip 


dky 

dky 

dky 


de 

777 

dp 


dkz 

dkz 

dk^ 

V 

de 

di} 

dp / 


(18) 


Note that J{e,'d,ip) = J ^{kx, ky, kz) with {e,'d,ip) expressed by Eq. (14). It is now definitely 
allowed to approximate the Dirac delta distribution by a Gaussian in analogy to Eqs. @ 


and (10) 


'BZ 


d^k ^(k) J(e(k)) 

de [d-d [ dp) \J{e,i^,p))\-g{e,i^,p))- 




exp — 




= / 9(k) 


TTcr 


exp 


e2(k) 




(19) 
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where in the last step the integration variables are changed back again to an integration 


over the wavevector using Eq. (14). This is exactly what we wanted to show in Eq. (13). 


One must keep in mind that the Jacobian determinant J{kx,ky,kz) given by Eq. (15) is 
zero for special points k = {k^, ky, kz) where Vke^(k) = 0. This is the case for the T-point and 
usually for points on the Brillouin zone boundarj^, and even the transformation (14) could 


be not continuously differentiable for special points. Then, the reverse transformation (16) 


used in Eq. (19) is not dehned anymore in a rigorous mathematical interpretation. However, 


these problems arise because of two idealizations, the long-time idealization and the inhnite- 
solid idealization, and the following remarks have to be considered: 

1. In a physical interpretation the Dirac delta distribution appearing in Fermi’s golden 
rule, Eq. ([^, is only a long-time idealization which should be replaced by a sin(x)/a;- 
function for realistic physical calculations. However, this would yield time-dependent 
rates which is usually not desired. 

2. For a numerical calculation an inhnite periodicity of the lattice is assumed (inhnite- 
solid idealization), k-points of this numerical calculation only sometimes coincide ex¬ 


actly with a point where the reverse transformation (16) is not dehned and the k-point 


grid can always be shifted so that there is no point where the reverse transformation 
is not dehned. For an arbitrary g{k) it is not clear that one gets a correct result when 
omitting these k-points. In a real solid in the ground state only a hnite number of en¬ 
ergy levels are occupied. These energy levels do not correspond to states with dehned 
wavevectors k. In the numerical treatment of these hnite systems the energy levels 
are approximated by the energies of a lattice with inhnite periodicity at a number of 
discrete points on a k-point grid. For sufficiently large systems the result must be 
independent of the detailed choice of the k-point grid. If we know the critical points, 
we can shift the k-point grid in such a way that it does not include the above dehned 
critical points (see also Ref. [IB]) . One could argue that a k-point very close to a critical 
k-point could yield an extremely large contribution because the Jacobian determinant 


(18) appearing in Eq. (17) may be very large for this point. This means 


that the second line of Eq. (19) would not be a good approximation for a choice of the 


k-point grid which contains points very close to critical points, but the near-equality 


of the hrst line with the third line of Eq. (19) still holds because for the transition 













between the first and the third line the integration variables have been changed back 
again and this corresponds to the multiplication with = J{kx,ky,kz), so 

altogether, a possibly large value of does not matter. However, it is ex¬ 

tremely complicated to identify all critical points and therefore it is not clear whether 
a chosen k-point grid contains critical points or not. In order to avoid these cum¬ 
bersome investigations one can also do the following: One performs calculations for 
denser and denser grids and/or for shifted and rotated grids and compares the results. 
If the results are very similar, this means that the grids either do not contain critical 
points or that the critical points do not make a big contribution so that they do not 
falsify the results, in agreement with Ref. [171 


III. PRACTICAL HINTS FOR THE APPROPRIATE CHOICE OF THE SMEAR¬ 
ING PARAMETER 


The appropriate choice of the smearing parameter a appearing in Eq. IT is crucial for 
the correct numerical calculation of the Boltzmann rate equation. The appropriate choice 
of a depends on two quantities: 


1. First, the smearing parameter a depends on the energy scale of the involved 
(quasi)particle which may be, e.g., a photon, phonon, magnon, plasmon (see Sec. |^. 
The energy scale for a phonon is in the order of some mRy (about 40 meV) and for a 
magnon the energy scale is a factor 10 larger than for a phonon. The energy scale for a 
plasmon is much larger, in the order of 700 mRy (about 10 eV). An appropriate choice 
of the smearing parameter is in the same order as the energy scale of the involved 
(quasi)particle. 

2. Second, the smearing parameter a depends on the grid spacing. A typical ansatz is 
a = p/Ni where Ni is the number of k-points in one direction and the total number 
of k-points is Nf. For a hxed proportionality constant p it is guaranteed that the 
smearing parameter is the smaller the larger Ni. 


In the following we discuss the choice of a for the case of electron-phonon scatterings. In 
Ref. |3j the smearing parameter is hxed to 15 meV (about 1 mRy) for the numerical calcu¬ 
lation of the electron-phonon Boltzmann rate equation. In Ref. |H we tested our numerical 
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FIG. 1: Rate of the magnetic moment change per atom dM/dt{ts) vs. number of k-points in 
one direction Ni for different smearing parameters a. For iron and an excitation temperature of 
Te = 2000 K. 

results of the electron-phonon Boltzmann rate equation for many different grids and smear¬ 
ing parameters. In this publication we considered the case of ultrafast demagnetization, see 
Sec. HI Among other quantities we calculated the rate of the magnetic moment change per 
atom dM/dt{ts) for a time tg (see Eq. (14) of Ref. Hj). tg is the time after the laser pulse 
irradiation where the electron system has thermalized, i.e., the electron distribution can be 
described by a Fermi-Dirac distribution with the electron temperature Tg. Fig. [T] shows the 
rate of the magnetic moment change per atom dM/dt{ts) for Fe and an electron temper¬ 
ature of Tg = 2000 K. We calculated dM/dtitg) for different grids (number of k-points in 
one direction A^i) and for different smearing parameters a. One can see that the results for 
dM/dt{tg) depend hardly on the chosen grid and on the chosen smearing parameter except 
for Ni = 10. Therefore, the above discussed critical points do not falsify our results and 
the smearing parameter is in the right order of magnitude (about several mRy). Of course 
it is trivial that increasing the number of k-points increases the convergence. By Fig. [^we 
just want to show that the results depend only very slightly on the specihc choice of a if the 
number of k-points is above a certain value. 
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IV. CONCLUSIONS 


Scattering processes of electrons at periodic perturbations are very important in solid 
state physics and also a correct numerical treatment is crucial for the quantitative analysis 
of scattering processes in many research activities, e.g., spintronics. In quantum mechanics 
Fermi’s golden rule is used which contains the Dirac delta distribution. The Dirac delta dis¬ 
tribution is usually replaced by a Gaussian in order to integrate numerically the Boltzmann 
rate equations on a hnite grid of k-points. It is not obvious from the very beginning that 
this numerical treatment is correct since the Dirac delta distribution is not a function and 
the smearing variable differs from the integration variable. We have shown in the present 
article that this procedure is in general correct. There are special k-points for which it is 
in principle not justihed to replace the Dirac delta distribution by a Gaussian, however we 
have given mathematical and physical arguments why this procedure is nevertheless a good 
approximation for the integration of the Boltzmann rate equation and should not falsify the 
results, at least for three dimensions. It is not clear whether the same holds also for d=2 
or even for d=l. In conclusion, the naive replacement of the Dirac delta distribution by a 
Gaussian gives in general correct results for the Boltzmann rate equation but this has to be 
checked for denser and/or for shifted and rotated grids in order to avoid wrong contributions 
from the above described special k-points where the replacement is critical. 
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